Confined polymer nematics: order and packing in a nematic drop 
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We investigate the tight packing of nematic polymers inside a confining hard sphere. We model 
the polymer via the continuum Frank elastic free energy augmented by a simple density dependent 
part as well as by taking proper care of the connectivity of the polymer chains when compared with 
simple nematics. The free energy ansatz is capable of describing an orient at ional ordering transition 
within the sample between an isotropic polymer solution and a polymer nematic phase. We solve 
the Euler-Lagrange equations numerically with the appropriate boundary conditions for the director 
and density field and investigate the orientation and density profile within a sphere. Two important 
parameters of the solution are the exact locations of the beginning and the end of the polymer chain. 
Pending on their spatial distribution and the actual size of the hard sphere enclosure we can get a 
plethora of various configurations of the chain exhibiting different defect geometry. 

PACS numbers: 64.60.an, 64.70.km, 64.70.pj, 64.70.mf, 64.70.Nd, 61.30.Dk, 61.30.Jf, 61.30.Pq, 61.30.Vx 



7d 

o 
o 



> 

p 

in 

o 
o 



't 



I. INTRODUCTION 

Packing of DNA within simple viruses has recently at- 
tracted a lot of attention from the physics community [l] 
since it appears that many if not all processes connected 
with the bacteriophage DNA injection are governed by 
biologically unspecific physical mechanisms. Cryomi- 
croscopy of simple viruses, such as bacteriophages T7 
[2], epsilonlS |3] and 029 [4], indicates that at elevated 
densities DNA appears to be wrapped into a coaxial in- 
verse spool, with pronounced ordering and high density 
close to the capsid wall that both appear to decay close to 
the center of the capsid. Addition of polyvalent count eri- 
ons such as the tetravalent spermine can induce a toroid 
formation inside the capsid reminiscent of the toroids ob- 
served in vitro [5 . These toroid-like packings seem to be 
observed in several bacteriophages but recent studies of 
Leforestier et al. [6l[7] on cryomicroscopy of T5 bacterio- 
phage indicate that more complicated packing geometries 
can also be realized such as nematic monodomains sep- 
arated by defect walls that in general do not conform to 
the inverse spool paradigm. Within this paradigm tight 
packing allows DNA to act like a coiled osmotic spring 
piled up against the inner surface of the capsid ready to 
release its osmotic and elastic energy through the portal 
complex on docking onto a bacterial wall |8]. 

The energetics of packaging and ejection of DNA in 
bacteriophages has been treated on various levels of ap- 
proximation [9Hl6] based on the Odijk - Gelbart inverse 
spool model in which the DNA chemical potential within 
the viral capsid is decomposed into an interaction part 
and a curvature part. This decomposition can describe 
the ejection of the genome reasonably when compared 
with osmotic stress experiments JTl [18] . Leaving aside 



simulations of the genome packing within the capsid 
[191 EQ] that rely on other sets of assumptions which we 
will not address in what follows, the different theoreti- 
cal approaches rely on the additivity ansatz for the total 
free energy of DNA packing confined to a spherical cap- 
sid. They assume that the free energy is composed of two 
parts: the curvature energy of DNA that is forced into the 
confines of the capsid, as well as DNA-DNA interaction 
free energy consistent with osmotic stress experiments in 
the bulk [21] |22] . In physical terms this additivity ansatz 
constitutes the basis of the inverse spool paradigm. 

The form of the elastic curvature energy is known, 
though some recent work might indicate additional sub- 
tleties that are usually not considered [23 . Its form, pro- 
portional to the square of local DNA curvature, follows 
from the standard Euler-Kirchhoff model of an elastic 
filament or in its continuum form from the Frank elastic 
energy. Though this model contains some subtle fea- 
tures due to the strong interhelical forces between the 
segments of the molecule [21 , it nevertheless appears to 
be a consistent description of DNA on mesoscopic scales 
[24 . The parameters of the Euler-Kirchhoffian model of 
DNA, such as its persistence length, are well established 
and have been measured by a variety of methods with 
satisfactory consensus among the results [25 . One can 
formulate the curvature elastic energy either on a single 
molecule level based on the Euler-Kirchhoff expression for 
the bending free energy or indeed on a continuous level 
with a volume distribution of polymer segments where 
the bending free energy stems from a general ansatz of 
the Frank elastic free energy [26] . 

The DNA-DNA interaction energy is less well known 
and its interpretation less well established. It can be 
measured directly in osmotic stress experiments 23 ^^^ 



can be deconvoluted into a longer ranged electrostatic 
contribution [21 and a shorter ranged hydration compo- 
nent [28]. Both of them have been quantified in terms 
of magnitudes and decay lengths [29l [30] . The variable 
directly deduced from experiments is thus the osmotic 
pressure in DNA arrays in the bulk and one can simply 
formulate the equation of thermodynamic equilibrium of 
encapsidated DNA [26] directly in terms of the measured 
osmotic pressure [30'-'32^ rather than in terms of theoret- 
ical polyelectrolyte models [9, 10, 33 or in terms of semi- 
empirical chemical potential expressions [12] [131 HSl US] • 
This approach yields a consistent DNA density profile 
within the capsid as well as the DNA encapsidation equa- 
tion of state, i.e. the dependence of the fraction of en- 
capsidated DNA on external osmotic pressure, that can 
be directly compared with experiments [TT] [18] . 

There are two major drawbacks and inconsistencies 
within this type of approaches. First of all these mod- 
els lead to a DNA-depleted region on the central axis of 
the capsid but do not take into account the possibility 
that within this depleted region the DNA is disordered 
and can not be described simply and solely in terms of 
its elastic free energy. On top of that one usually does 
not solve directly for the polymer director profile within 
the capsid enclosure but assumes an inverse spool ansatz 
form from the start. For a completely consistent descrip- 
tion one needs a free energy ansatz that would allow for a 
nematic-isotropic transition of the DNA solution as well 
as describe the director field configuration consistent with 
the distribution of the free ends of DNA within the cap- 
sid. 

In what follows we will thus venture to set up a model 
of polymer nematic ordering within confined enclosures 
that would take into account the coupling between the 
density field and the director field of a polymer nematic 
as well as the possible role of the ends of the polymer 
chain that act as nucleators of defects in the director 
field. 

We start our analysis with a nematic ordering free en- 
ergy appropriate for the case of long semifiexible poly- 
mers and derive the free energy contributions that are 
due to the well known coupling between the polymer den- 
sity and director fields. This allows us to calculate the 
equilibrium director and density profiles within a spher- 
ical enclosure without any additional assumptions. We 
then investigate various packing configurations of con- 
fined polymer nematogens in order to asses the value of 
the inverse spool paradigm. 

Though our calculations are motivated by the prob- 
lem of DNA packing within viral capsids, our results 
are just as relevant for a completely general case of 
confined semifiexible polymer ordering within spherical 
enclosures. While identification of our approach with 
the physics of capsid- confined DNA is only approxi- 
mate since the persistence length of DNA is on the or- 
der of the capsid size, the appropriateness of the model 
would become more accurate for larger radii of the enclo- 
sure. It also addresses a previously seldom studied and 



thus poorly understood fundamental problem of confined 
polymer nematics. In any case, we are convinced that our 
mesoscopic approach nicely complements molecular sim- 
ulations and molecular mechanics that has been used in 
the context of DNA packing within viral capsids before. 



II. THEORETICAL DESCRIPTION 

Since DNA is a long polyelectroyte molecule its con- 
nectivity introduces additional features into a consistent 
continuum description of its nematic ordering that are 
not addressed by the standard approach of the liquid 
crystal physics [34 . Based on previous work by Kamien 
et al. [35] we thus first introduce and formulate the con- 
cept of the polymer current and then construct an ap- 
propriate nematic ordering free energy that we solve in 
confined geometry of a spherical capsid. 



A. Polymer "current density" 

In a polymer nematic liquid crystal, splay deforma- 
tion becomes progressively more expensive with increas- 
ing chain length, as it imposes costly local changes in 
polymer density [35, 36 . In the continuous limit of long 
chains this coupling between the polymer density and 
orientational fields is described by an analogue of the 
continuity equation for the nematic director field n [37], 



V . (psu) = 0, 



(1) 



where ps is the surface number density of chains crossing 
the plane perpendicular to the director field [38j. Eq. (!]) 
can be interpreted as the continuity condition for a poly- 
mer "current density" j = pgii. The only difference be- 
tween Eq. ^ and the usual continuity equation is that 
in this case j does not describe a rate (there is no time 
derivative involved in its definition), i.e., we are observ- 
ing the number of chains perforating the perpendicular 
plane rather than the number of particles crossing it per 
unit time. For the same reason the time-dependent term 
(■^ or similar, where p is the volume density) is absent. 
The analogy comes fully into life if one relaxes the con- 
dition |n| = 1 and takes into account the variable degree 
of nematic order. Let us stress that at this stage we will 
consider polar ordering rather than quadrupolar nematic 
ordering. It is not yet understood how to construct a re- 
placement for the polymer current in case of the nematic 
order tensor. In principle, there is no polar ordering in 
a regular nematic. Yet in the case of the polymer ne- 
matic its definition is still useful. Here it can be defined 
because of the connectivity of individual long chains, in 
opposition to the case of shorter nematogens. For suffi- 
ciently short times (such that the diffusion of ends is not 
effective) or if the ends are pinned, the ordering is polar, 
if we identify the beginning and the end of the chain. 



The description of nematic ordering will be introduced 
through a non-unit nematic order vector a, 



a = {cos 6) n. 



(2) 



It is defined in a given hydrodynamic volume V by the 
expression 



1 




^ ^(y,( [dhcosOi])n, (3) 



where the integrals go over the length of the ith. chain 
within the hydrodynamic volume and the sum goes over 
all the polymer chains whose parts can be found within 
the hydrodynamic volume; 6i is the angle between the 
local tangent on a chain and the average vicinal direc- 
tor field represented by n, and () is the thermodynamic 
average. 

Let us now calculate the "fiux" of chains through a 
surface perpendicular to the average vicinal director field 
n, i.e., the number of chains perforating this surface per 
unit area. On average, a subunit i of the chain perforates 
the surface if the distance between its center and the 
surface is smaller than ^ J^ dli cos Oi on either side of the 
surface; io is the subunit length. Averaging over all chain 
configurations in the hydrodynamic volume we obtain 

(EJ,d^.cos^.) {EJudk) 

Here we introduced the volume number density of sub- 
units p and used the definitions Q and ([3|, with (TV) = 
(}Zi 1) ^h^ average number of subunits in the hydrody- 
namic volume, (^i Ji.dk) /{N) is obviously the length 
of the subunit, Iq. Note that the length of the subunit 
may be chosen arbitrarily, in relation with the definition 
of the density, i.e., the product pio is independent of this 
choice. 

By definition this "fiux" density of chains through the 
surface is given by the normal component of the "current 
density", jn = j • n. We can thus define the polymer 
current density at any point within the sample as 



p^oa. 



(5) 



In the infinite chain limit where there exist no chain be- 
ginnings or ends, the number of chains entering a closed 
surface must be equal to the number of chains exiting it, 
hence 



V- 



0. 



(6) 



Moreover, if the chain beginnings and ends do exist, then 
the fiux through a closed surface is nonzero: it is positive 



if there are more beginnings than ends within the closed 
surface, and negative in the opposite case. Thus, the net 
density of the beginnings and ends of the chains acts as 
a source term and we can write down a complete local 
continuity requirement. 



Vj = p^ 



(7) 



where p^ is the volume number density of the beginnings 
{p^ > 0) and the ends {p^ < 0) of the chains. Eq. ^ 
presents a generalization of Eq. (IT]), taking into account 
variable ordering and source terms due to chain ends. 
It comes in the form of the usual continuity equation, 
except that, by construction, it does not contain the time- 
dependent term. 

The microscopic definition of the polymer current den- 
sity ([5| and the continuity equation ^ can be obtained 
as follows (for details see e.g. [39]). Let us first show that 
the vector 



t(r) 




dk rili)5^ir - r{k)) 



(8) 



is just the polymer current density (|5|. Above the sum 
goes over all the polymer molecules in the system, and the 
integrals go over the full length of each polymer; r(/^) = 
^j:^ is the local unit tangent vector of the chain. An 
insight into the nature of t is obtained by smoothing it 
out by integration over a hydrodynamic volume V as 

where the last sum and integral go over the chains and 
the segments within the hydrodynamic volume F, respec- 
tively. Furthermore, by writing the microscopic subunit 
density of the polymer chains in the form 



p{r) = UY.J^^dk5(^\r-r{k)) 



(10) 



where €o is again the length of the subunit, we obtain that 
the smoothed hydrodynamic subunit number density is 



1 
V 



d^r p{r) 



/dV /^ /"d^i,5(3)(r 



1 1 

VTojy 
1 1 

V'^0— \Ji, 



r{k)) 



E( / d^,: 



Equations (11) and ^ then give back exactly the defi- 
nition of the polymer current density ([5| with the iden- 
tification that the vector t, Eq. ([8|, is exactly the micro- 
scopic version of the current density j of the polymer. 

This identification of the polymer current density leads 
directly to the continuity requirement ^ that it needs 



(9) 



(11) 



to satisfy. This can be seen by evaluating the divergence 
of the microscopic current density vector, obtaining 



= - ((y2sHr-n{L))\ - (rsHr-rM))] ■ 

(12) 

Here we took into account that V — > — V^, noting that 
the argument of the Dirac delta function is r — r(/^). 
Also, r{li) • Vi = ^. The above result can be written 
alternatively as 



where 



V-j(r) = p+(r)-p-(r), (13) 



o+{r) = {Y,5'{r-rm) 



and 



p-ir) = {^5Hr-r,{L)) 



are the total densities of the beginnings and ends of the 



chains in the sample. Identifying p+ — p' 
thus formally recovered Eq. ([7|. 

This constraint, Eqs. (fTl) and (fl3|, on the director field 
in polymer nematics has been discovered by de Gennes 
and Meyer [37] and bears some similarity with the differ- 
ential form of the Gauss theorem in electrostatics. The 
divergence of the field is determined by the density of the 
sources. Here the orientational field j plays the role of 
the electrostatic field and the density of beginnings and 
ends of the chain play the role of positive and negative 
charges. 

In the limit of fixed nematic order and fixed density, 
only if there exists a proper mismatch of p^(r) and p~ (r) 
can there be a splay deformation of the polymer nematic. 
Eq. (13) presents a generalization of Eq. (IT]), taking into 
account variable ordering and source terms due to finite 
length of the polymer chains. 



B. Free energy 

To determine the equilibrium configuration of the di- 
rector and density fields we set up a free energy density 
following the Landau approach [34 . The appropriate 
variables in the polymer case are the complete director 
field a, describing the orientation and the degree of order, 
and the polymer density field p. Both fields are coupled 



by the continuity requirement ^. The conservation of 
polymer mass will be satisfied globally. The phase tran- 
sition will be controlled by the density (concentration) 
of the polymer as is the case for lyotropic nematic liquid 
crystals. 

Let us stress that for computational reasons all equa- 
tions must remain regular also for vanishing order, i.e., 
they must be expressed by the full vector a. Decomposi- 
tion of the form a = an, where a is the degree of order, 
would result in a singularity of the form times oo taking 
place in centers of defects, where Vn diverges while the 
degree of order vanishes. In contrast, a and its deriva- 
tives remain regular everywhere. Taking into account the 
definition ([5|, Eq. ^ is already of the correct form. 

In the elastic free energy, instead of using the usual 
Frank terms for splay, twist, and bend of the director. 



a new set of elastic terms must be used [40l [41] : 



(14) 



+ -L'3aiaj{diak){djak) + 



(15) 



L4^{eijkakdiaj) . 



Unlike the Frank elastic parameters Ki, the elastic con- 
stants L[ do not depend on the degree of order. To keep 
the number of elastic parameters at minimum, among 
all possible terms quadratic in the derivative we have re- 
tained only those non-vanishing in the limit of a fixed 
degree of order. Comparison of Eqs. ( 14) and ( 15 ) in this 



p , we have limit relates the Frank constants Kj to the constants L,-: 



Ki 


= a^L[+a^L'2, 


K2 


= a^L'i+a'^L'^, 


K3 


= a^L'i+a'^L's. 



(16) 
(17) 
(18) 

One observes that the dependence on the degree of order 
is different than in the case of the nematic tensor or- 
der parameter [t40J |41J , where the corrections that break 
the degeneracy of splay and bend are cubic in the de- 
gree of order. Note that all three deformation modes 
(splay, twist, and bend) are degenerate in the limit of 
small degree of order, as indicated by the leading order 
dependence Ki (x o? . This would not be the case if one 
just used a in the Frank expression (14). Yet it must 



be so, because the direction n necessary to distinguish 
between splay, twist, and bend, is not defined as a ^ 0. 
Generally, the Frank expression (14) can be used only 
perturbatively, i.e., for |a| ^ 1, whereas the expression 
(15) or similar can be used to describe the full range 



<a < 1. 

To establish a one-to-one correspondence between the 
two sets of elastic parameters, we make a further simpli- 
fication of our model free energy expression (15), while 



retaining full elastic anisotropy known to be significant in 
lyotropic liquid crystals. One of the possibilities is omit- 
ting the I/4 term. Hence, the minimal free energy model 



in this case reads 



C. Euler - Lagrange equations and their solution 






'-a- 
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-pCa^ 



pLi{diajf 



(19) 
(20) 
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-pL2{diaif + -pL^aiaj{diak){djak) (21) 
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n2 



di{pai) 



io 



X{P- Pof 
Lp{dipf, 



(22) 

(23) 

(24) 



where C is a positive Landau constant describing the 
isotropic-nematic phase transition and p* is the transi- 
tion density, 



Li = K2/ipoa^), 

L2 = {K^-K2)/{poa^), 

Ls = {K3-K2)/{poa''), 



(25) 
(26) 

(27) 



po is the bulk equilibrium density, x ^^^ ^p ^^^ the den- 
sity compressibility and the density variation correlation 
length. The nonlinear density factor in the first term of 
( 19 ) guarantees that the bulk nematic ordering stays lim- 
ited to |a| < 1. In the terms (19)-(21) of the total free 



energy, corresponding to the ordering and elastic parts of 
the free energy, we have taken into account the fact that 
they need to be proportional to the number of molecules, 
i.e., to the local density. 

The continuity requirement ^ is taken into account 
by means of the penalty potential ( |22| ) proportional to 
a coupling constant G (the arbitrary length £o has been 
absorbed in G). For most of the time, the density of 
beginnings and ends of chains, p^, will not be considered 
as a variable but as a fixed external parameter. 

What physical meaning can be attributed to the pa- 
rameter G that defines the rigidity of the continuity con- 
straint ^ , as this parameter does not exist if the continu- 
ity constraint ^ is satisfied exactly? One notes that for 
the case p^ = 0, there emerges a splay relaxation length 
scale Ip = \JG Ix that controls the rigidity of the conti- 
nuity constraint. Splay deformation of the director field 
at length scales much larger than Ip does not bring about 
any substantial density variation, i.e., at this scale splay 
and density are decoupled. On the other hand, on length 
scales much shorter than /p, the coupling between splay 
and density becomes strong. This makes sense physi- 
cally as on large length scales the (weak) divergence of 
the polymer current is compensated by a spontaneous 
rearrangement of chain ends - which we are not taking 
into account - while this is not effective at shorter length 
scales. It seems natural that Ip should be in direct rela- 
tion with the chain length [35] , 



We will find the equilibrium configuration of both con- 
stitutive fields, i.e. the density field p and the non-unit 
nematic director field a, by minimizing the free energy at 
the constraint of global polymer mass conservation. The 
corresponding functional is thus 



dF(/-Ap), 



with the constraint 



dV p = mo = const. 



(28) 



(29) 



where A is a constant Lagrange multiplier. The mini- 
mization will be performed by following a quasi-dynamic 
evolution of the director and density fields of the form 



daj 
dt 
dp 
di 



7- 



7 



= d, 



= d, 



df 



d{djai 
df 



_dl 

dai 

df 



+ A, 



(30) 



(31) 



d{djp) dp 

where 7 is a formal parameter defining the time scale. 
Eq. (31 ) shows that satisfying the constraint of mass con- 



servation is easy: at every time step one subtracts from 
the density the homogeneous field 

Ap=^ JdVp-mo (32) 

SO that the corrected density p — Ap satisfies the mass 
conservation (29). 



We use a tangentially degenerate boundary condition 
for the director, i.e., the director is everywhere parallel 
to the surface of the sphere, while it is allowed to rotate 
freely in the tangential plane. For the density we use free 
boundary condition, i.e., the normal (radial) component 
of the density gradient is zero. 

The initial condition is p(r) = po for the density field 
and a(r) = plus a small random perturbation for the 
director field. The equations are solved by an open 
source finite volume solver (OpenFOAM accessible at 
|http://www. openfoam.c om/ ) on a 50 x 50 x 50 cubic 
mesh (size of the box containing the sphere), which is 
deformed and refined near the surface of the sphere to 
define a smooth boundary. 



III. RESULTS 

In what follows we present the steady state solutions 
of the quasi-dynamic evolution of the director and den- 
sity fields corresponding to direct solutions of the Euler- 
Lagrange equation. The steady-state solutions obviously 
satisfy 

df dl 

dai 
df 







^'d{djai) 



(33) 







di 



df 
didjp) dp 



+ A, 



(34) 
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FIG. 1: The director and density fields in the case of no 
sources, p =0. The radius of the confining sphere is (a) 
10 and (b) 32 in units of the nematic correlation length. The 
average free energy density / of the configurations is indi- 
cated; it should be compared for systems of the same size. 
The director field is shown by random tracer paths, on the 
iso-contour the value of the nematic order is half the bulk 
equilibrium value. The density profile is shown in the cutting 
plane. 



with the appropriate mass constraint (29). Unless stated 



otherwise, in the numerical simulation we take the follow- 
ing dimensionless values for the elastic constants entering 
our model: A = C = 1, Li = 1, L2 = L3 = 0, G = 1, 
X = 1, Lp = 1. This means that in these units the 
nematic correlation length (the characteristic size of the 
defect core) is 1, as is the characteristic length Ip of the 
director-density coupling. The bulk equilibrium density 
is set to po = 1-5 and the nematic transition threshold 
density to p* = 0.5. The only remainin g pa rameter, -^0, 
merely scales the source density p^, Eq. (22), and can be 
set to ^0 = 1 without loss of generality. 

The tangential boundary condition for the director im- 
plies a disclination of the total strength +2 on the sur- 
face of the sphere [3T . Most usually this means a pair of 
disclinations with strength +1 each. Quite generally, due 
to the expensive splay deformation in nematic polymers, 
circular disclinations are preferred to the radial ones 
which would be favoured in a regular nematic because 
they are compatible with a cheaper three-dimensional 
configuration. 



A. No sources 

We first of all find the solutions for the case of no 
sources, p^ = 0, in a spherical hard enclosure. Obvi- 
ously, see Fig. ma), the geometry of the director field is 
pronouncedly toroidal of an inverse spool type (note that 
the director field is represented by trajectories of random 
tracers), while the density field shows a moderate deple- 
tion on the cylindrical axis. Most notably, the density 
depleted axial core is only partially orientationally disor- 
dered, i.e. the disclination line is avoided by the director 
spiraling becoming increasingly vertical in the core. The 
iso-contours of the nematic ordering (|a|) near the poles 



FIG. 2: Dipolar sources: p^ — I'm the green region, p^ — —1 
in the red region. The radius of the confining sphere is (a) 10 
and (b) 32. The key is explained in Fig.fl] 



show that the defects are predominantly point-like but 
not well isolated due to the small system size. 

The equilibrium profiles of the density and director 
field depend crucially on the size of the confining sphere. 
Making the sphere bigger, see Fig. [lib), leads to a pro- 
nounced localization of the axial density depletion range 
to polar regions, while the rest of the spherical volume is 
almost uniformly packed with the polymer as the director 
configuration is more spiral and less toroidal if compared 
to the case of tighter confinement. This can only come 
about at the cost of a deformed director field that devel- 
ops two opposite vortices (with strength +1 each) in the 
polar regions. 



B. Fixed sources 

Introduction of sources further modifies the constitu- 
tive fields. We first of all take two oppositely "charged" 
localized sources, p^ = 1 and p^ = — 1 located sym- 
metrically at the two poles, see Fig. [2] Since the sources 
correspond to oppositely " charged" ends of the chain we 
can refer to this configuration as a dipolar configuration. 
The introduction of the two sources has only a limited 
effect on the configuration, as the configurations in Fig.[l] 
already exhibit a cylindrical symmetry (with the distinc- 
tion that there the direction of the symmetry axis is se- 
lected spontaneously). The major difference is that now 
the two surface defects become more radial, and hence 
the configuration is less toroidal. This is expected since 
now the sources allow more splay around the defects. The 
pulse-like behaviour of the density around the point de- 
fects is a signature of the step-like profile of the source 
densities. 

For completeness we also show, see Fig. [3) the case of a 
configuration characterized by a single localized source, 
while a source density of opposite sign is evenly dis- 
tributed over the whole system such that the total source 
is zero. This would correspond, e.g., to the polymer chain 
with one end fixed to the wall and the other end free -a 
situation often encountered while DNA is being packed 
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FIG. 3: Single source, p = 1 in the green region, a corre- 
sponding negative p is evenly distributed elsewhere so that 
the total source is zero. The radius of the confining sphere is 
(a) 10 and (b) 32. The key is explained in Fig.fl] 



FIG. 5: Density variation cheaper, no sources. The radius of 
the confining sphere is (a) 10 and (b) 32. The key is explained 
in Fig.[l] 





(a)/: 



-0.06910 



(b)/: 



-0.08807 



FIG. 4: A p =1 (green) and p = — 1 (red) source placed 
on the z and x axes. The radius of the confining sphere is (a) 
10 and (b) 32. The key is explained in Fig. IT] 



into a capsid. While the single source configurations are 
quite similar to the dipolar case, it is instructive to ob- 
serve the difference between the two defects of the single 
source configurations. The defect located at the end op- 
posite to the source is almost entirely circular, in contrast 
to the defect located at the source. 

An oblique placement of the localized p^ = 1 and 
p^ = —1 sources, Fig.jij shows a particularly significant 
dependence on the size of the confining sphere. In the 
case of strong confinement, corresponding to a smaller 
radius of the sphere, the configuration is almost unaf- 
fected by the presence of the sources. While in the case 
of weaker confinement, corresponding to a larger radius of 
the sphere, the director configuration is completely gov- 
erned by the two sources and shows the same symmetry 
as the positions of the sources. 



C. Compressibility and density correlation length 
effects 



The director configurations and packing density 
change significantly if we make density changes, as quan- 
tified by the compressibility modulus x ^^d the density 



correlation length Lp, less costly in energy terms. To 
demonstrate this we choose rather drastically x = 0-1 
and Lp = 0.33, which nevertheless still yields p(r) > ev- 
erywhere within the enclosed volume. The severely con- 
fined configuration with no sources, Fig.jsja), is in this 
case purely toroidal and exhibits a clear axial disclina- 
tion line, while the density field shows a strong depletion 
around the axis of cylindrical symmetry. Moreover, in the 
case of a bigger confining sphere, Fig.jsjb), the polymer 
is completely depleted from a part of the available space, 
breaking the polar symmetry. It appears that stronger 
confinement operating for smaller sphere radii induce a 
toroidal packing irrespective of all the other terms in the 
free energy. This finding should be particularly relevant 
for the DNA packing within viral capsids and points to 
a universality of the inverse spool configuration. 

The configurations with the dipolar sources show a pro- 
nounced central density depletion region and a stronger 
nematic order reduction in case of severe confinement cor- 
responding to smaller sphere radii, Fig.JGJ^a). The case of 
a bigger sphere radii again leads to purely polar defects, 
Fig.lHFb), with a density variation showing a local spher- 
ical symmetry. In regular nematics, the latter director 
configuration is readily encountered and is known as the 
bipolar director configuration [42^. 

The single source cases, Fig.[7| again show a strong ax- 
ial polymer depletion and decrease of ordering (a), as well 
as almost complete depletion around a single pole within 
the sphere (b). The configurations with oblique sources, 
Fig. [S] show similar features as before, plus a complete de- 
pletion in the part of the bigger sphere which is somewhat 
different to Figs. |5] and [7| In the latter cases the depletion 
conveniently removed one of the defects, whereas in the 
present case it is the most sluggish part of the director 
field which is depleted. 



D. Multiple sources 

We have analyzed the case of no sources, the case 
of symmetrically and asymmetrically placed dipolar 







(a)/: 



-0.07308 



(b)/: 



-0.08925 



(a)/: 



-0.08740 



(b)/: 



-0.08727 



FIG. 6: Density variation cheaper, dipolar sources as in Fig.[2] 
The radius of the confining sphere is (a) 10 and (b) 32. The 
key is explained in Fig.[l] 





(a)/: 



-0.07458 



(b)/: 



-0.09182 



FIG. 7: Density variation cheaper, single source as in Fig.[3] 
The radius of the confining sphere is (a) 10 and (b) 32. The 
key is explained in Fig.[l] 



sources, and the case of a single source. We now deal 
with less symmetric source configurations corresponding 
to multiple sources. As an example we show, Fig. [9) a pair 
of dipolar sources in parallel and antiparallel arrange- 
ments in the case of less pronounced confinement within 
a bigger sphere. The parallel configuration, Fig.^a), 
is naturally formed by two halves with two depletion re- 
gions for space filling reasons. The antiparallel configura- 





(a)/: 



-0.07374 



(b)/: 



-0.08988 



FIG. 8: Density variation cheaper, sources placed on the z 
and X axes as in Fig. [4] The radius of the confining sphere is 
(a) 10 and (b) 32. The key is explained in Fig.[l] 



FIG. 9: A pair of dipolar sources, p = ±1, (a) parallel and 
(b) anti-parallel. The radius of the confining sphere is 32. 
Note in (b) that the director field inside the closest green 
region is defect-free, yet containing a well visible splay defor- 
mation. The key is explained in Fig.^ 




(a) 



(b) / = -0.08967 



FIG. 10: A pair of antiparallel source dipoles, Fig.lMb): (a) 
early stages of the evolution showing a more graspable symme- 
try than the actual equilibrium state in Fig.jojb). In (b) the 
equilibrium configuration in the presence of elastic anisotropy 
as discussed in Sec. lIII El is shown. Note that in this case all 
four source regions are populated by defects, in contrast to 
Fig.l9|b). The radius of the confining sphere is 32. The key 
is explained in Fig.[l] 



tion, Fig.|9|b), is less regular and seems to be quite frus- 
trated, yet the average free energy density is only slightly 
larger than in the parallel case. In the early stages of the 
evolution, the antiparallel configuration appears to be 
more tractable. Fig. 10 'a), but is then transformed into 
the equilibrium state of Fig.loFb). 



E. Effects of elastic anisotropy 

In the calculations presented above we have assumed 
everywhere that L2 = Ls = 0, i.e. we assumed a one 
constant approximation of the nematic elasticity. This 
limitation should not be too severe since splay is inher- 
ently different due to the coupling to the density, but one 
nevertheless needs to check whether there are any sub- 
stantial variations in the equilibrium configurations if the 
elasticity is anisotropic, in particular in those exhibiting 
the depletion regions. 




where x± is a source density compressibility. It is re- 
quired that 



(a) / = -0.09390 



(b) / = -0.09219 



FIG. 11: Anisotropic Frank elasticity, Li — L2 — 0.5, L3 = 1: 
(a) no sources, (b) oblique sources. The radius of the confining 
sphere is 32. The configurations are to be compared to those 
on Figs.pTb) andlSTb). The key is explained in Fig. IT] 



We thus assume now that the elastic anisotropy is 
present and we fix the elastic constants to Li = L2 = 0.5, 
1/3 = 1. Obviously, Figs. 11 andllOfb), the general fea- 
tures of the equilibrium nematic order and density con- 
figurations presented so far, in particular the nontrivial 
cases with the depletion, are robust with respect to the 
introduction of elastic anisotropy. 



In Fig. 11 the con- 
figurations with (a) no sources and (b) oblique sources 
are shown while the other parameters are the same as in 
Figs.[5]and[4| 

Most notably. Fig. 



11 ^a), the polymer is again com- 



pletely depleted from a part of the available space. There 
is no strict inverse-spool ordering like in Fig.jsjb), i.e., 
the axial defect line is substituted for the point defect 
via spiraling. Such elastic-anisotropy-dependent struc- 
tural transitions are well known in regular nematics [42 . 
In the case of oblique sources, Fig.jTTJb), the director 
field is much more irregular than that of Fig.lsFb), but 
the main feature - the depletion - is still there. One can 
thus conclude that with anisotropic elasticity only details 
of the director field change, while the density depletions 
and the overall geometry of the nematic order remain the 
same as in the isotropic case. 



F. Self-consistent distribution of sources 

So far we have studied cases without or with a fixed dis- 
tribution of sources. Let us finally check whether the dis- 
tribution of sources can be determined self-consistently 
together with the director and density fields. In other 
words, p^{r) will be treated as an additional variable 
rather than a fixed field as before. Of course, an addi- 
tional penalty term for source variations must be added 
to the free energy density (19)-(24). The simplest possi- 



bility is a quadratic penalty potential 



fs = 7:X±{P ) , 



dVp^ =0, 



(36) 



therefore an additional Lagrange multiplier is introduced. 
The procedure of satisfying this constraint is fully anal- 
ogous to satisfying the mass constraint, Eqs. (28), (29), 
(IsTl, and (l32|). 



(35) 



In the trivial case where x± is zero, the density of 
sources simply adjusts such that it exactly compensates 
for any divergence of pa and thus the director-density 
coupling energy (22) is zero identically, i.e., this corre- 
sponds to putting G = 0. For any nonzero Xi, however, 
the source density is penalized and cannot adjust freely. 
The configurations obtained with x± = 0-5 are shown 
in Fig.[T2j The source density is well-defined and is typi- 
cally dipolar, favouring the bipolar director configuration 
in the small sphere. Fig. 12 ^a), over the axially depleted 
one in Fig.lsja). These two configurations can be com- 
pared directly, as all parameters, apart from the sources 
variation, are the same. In the bigger sphere, Fig.[l2|b), 
one half of the bipolar director field is traded for the de- 
pleted region which absorbes the two point defects of the 
bipolar director configuration. The director field resem- 
bles the one in Fig.jsj^b), but there the defects appear 
to be forced by the fixed localized sources to stay away 
from the depletion. One can verify that in Fig. 12 a) and 
(b) the average free energy density is indeed lower than 
in any other configuration that can be compared directly 
with it. 

If, however, the source penalty gets larger, x± = 1 in 
Fig.[T3J then the source density falls towards zero rapidly 
and more so for the tighter confinement. There appears 
to exist a structural threshold, above which the source 
density is zero exactly. The verification and quantifica- 
tion of this threshold is not within our aim at present. 
The configurations in Figs.[l3]^a) andjsja) are identical, 
as is the average free energy density. Moreover, the fact 
that the average free energy density is slightly (!) higher 
in Fig.fl3[a) tells us that the configuration is not yet com- 
pletely relaxed to equilibrium and so the source density 
will fall further. On the other side, the average free en- 
ergy in Fig. 13 ^b) is lower than in Fig.lsjb), indicating 
that the source density is actually nonzero as presented 
in the figure. Apart from that, the configuration is very 
similar to Fig.lsjb). 



IV. CONCLUSIONS AND DISCUSSION 

Ordered structures containing a small number of di- 
rector defects are typically encountered in regular ne- 
matics when subject to a tight confinement that exhibits 
a sufficiently strong director anchoring at the confining 
walls [42 . Configurations like those in Figs.lsja) and|6|b) 
are typical for regular nematics in spherical confinement 
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(a) / = -0.07727 



(b) / = -0.09318 



FIG. 12: a self-consistent source field: p is treated as a 
variable; x± — 0-5, density variation cheaper. The radius of 
the confining sphere is (a) 10 and (b) 32. The key is explained 
in Fig.^ The source field is shown in the additional red-green 
cross section. 





(a)/: 



-0.07630 



(b)/: 



-0.09288 



FIG. 13: a self-consistent source field: p is treated as a 
variable; x± — ^^ density variation cheaper. The radius of the 
confining sphere is (a) 10 and (b) 32. The key is explained in 
Fig.[l] The source field is shown in the additional red-green 
cross section. 



(e.g. nematic droplets). The variable density and its cou- 
pling to the splay deformation in case of nematic poly- 
mers opens the door to a variety of new features and 
configurations. The spontaneous appearance of depleted 
regions (phase separation) is particularly interesting in 
this context. On the other hand we must not forget that 
by using the polar ordering instead of the usual nematic 
quadrupolar ordering, we give up a class of configurations 
containing defects of half-integer strengths [43] . 

T5 cryoelectron microscopy experiments [5^ show that 
in case of a loose packing inside a bacteriophage cap- 
sid, the DNA toroids condensed with spermine {sp^^) 
often occupy only a part of the available volume rather 
than being spread over the entire capsid. Contrary to 
the toroids observed in the bulk [4¥, the toroids con- 
fined to a viral capsid often show no polar symmetry. 
The nematic and density field configurations obtained in 
our work, containing pronounced depletion regions, show 
that the partial packing scenario could be in principle at- 
tributed already to a macroscopic free energy of the type 
(19)- (24), as an alternative to the strictly cohesive energy 



gies are inferred from experiments on the bulk polyvalent- 
counterion condensed DNA j 45l - [47] where the nematic or- 
der as well as the DNA density are homogeneous within 
the sample. In viruses the DNA-DNA cohesion mediated 
by polyvalent counterions such as sp^^ O [6], coupled 
with strong spatial confinement and nematic elasticity, 
might play an altogether different role leading to struc- 
tures diverging from a simple inverse spool paradigm with 
a pronounced lack of polar symmetry. Our numerical re- 
sults certainly point to the same conclusion: we obtain 
depleted regions without polar symmetry in the cases 
with larger confining spheres, Figs.[5J [Tj [8J The same 
conclusion regarding the absence of polar symmetry in 
the encapsidated toroidal aggregate holds also in the case 



of anisotropic elasticity. Fig. 11 , where the polymer is ex- 
pelled from one polar region but not from the other. We 
note here that our simplified macroscopic free energy, 
(19)- (24), can not capture all the subtleties of the inter- 



segment interactions of a real nematic polymer such as 
DNA [211 [22]. 

Recent elucidation of the organization of the DNA 
packing in the T5 capsid at various DNA densities as 
it undergoes a series of phase transitions on decreasing 
the DNA density [7 attests to the possible minor role of 
inverse spool states in viral DNA packing [48]. Though 
our results definitely point to the existence of states with 
partial depletion regions in polymer packing with no po- 
lar symmetry, they also indicate that the inverse spool 
family of states is nevertheless quite robust within our 
model assumptions. This finding is completely consis- 
tent also with recent molecular simulations, see below 
^49] . Strong confinement for smaller confining spheres 
appears to induce toroidal packing irrespective of all the 
other terms in the free energy and points to a univer- 
sality of the inverse spool configuration. The observed 
non-inverse spool states might thus originate in other 
types of interaction energy e.g. between the DNA and 
various molecular moieties or charged groups distributed 
along the inner surface of the capsid, effects which are 
not considered in this work. 

Experiments on DNA ejection from partially filled cap- 
sids also shows a profusion of monodomains of nematic 
order separated by dislocation, twist and bend walls [7]. 
Though we have explored in relative detail the parame- 
ter space of our numerical solutions of confined polymer 
nematic packing we were unable to detect packings with 
domain structure of the type observed in these experi- 
ments. The same is true for the experimentally observed 
sequence of liquid crystalline phases on ejection of DNA, 
starting with hexagonal and ending in isotropic [7]. Our 
free energy and nematic order description is of course in- 
sufficient for the description of these effects. A serious 
limitation of our model might be the simplified form of 
the density dependence of the free energy which is basi- 
cally reduced to the quadratic deviation term, Eq. ( 23 ) . 



between DNA segments. Cohesive interaction free ener- 



The osmotic pressure measured directly for bulk DNA 
[21] [22] certainly shows the inadequacy of this simplifica- 
tion and would in fact lead to a very complicated density 
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variation of the free energy with pronounced salt and 
polyvalent counterion concentration dependence. Unfor- 
tunately the details of this dependence are not known 
for a wide range of parameters and are understood even 
less on a molecular level. Our choice for a simplified 
version of the density dependence is thus in this respect 
fully vindicated by the blind spots in our understanding 
of the fundamental interactions between DNA molecules. 
It seems appropriate that theoretical formulations of the 
DNA viral packing problem remain as close to the exper- 
imentally determined variables as possible [26 . 

Our numerical solution of the confined polymer ne- 
matic model also provides an interesting alternative to 
simulations of DNA packing within icosahedral viruses 
j49] . Since the seminal simulation work of Kindt et al. 
[TT] inverse spool-like structures have been obtained on 
various levels of sophistication. Typically in a simula- 
tion a semifiexible chain, possibly with chiral twist in- 
teractions as in the work of Spakowitz and Wang [20] , is 
endowed with different types of self-interaction that con- 
tains repulsive screened electrostatic and/or attractive 
ion-correlation or hydration contributions. The DNA 
packaging starts at the inner capsid surface, typically 
with a helical winding along the inner surface of the 
sphere-like capsid [50, 51 , and then proceeds inward to- 
wards a more disordered core using all available inner 
space [I9l [20l |52] . Sometimes the capsid is modeled with 
more detail than in the spherical model, including also 
icosahedral symmetry of the viral capsids [53^. Details of 
the interaction potential also seem to be very important. 
Presence of intersegment attractions in DNA-DNA inter- 
actions has been shown to lead to configurations that do 
not conform to the inverse-spool paradigm. Forrey and 
Muthukumar were the first to observe a folded toroid 
state in DNA packing with attractive self-interactions 
within an icosahedral enclosure [54]. The same pack- 
ing configuration was later seen also in [55]. The spools 
and folded toroids observed in these simulations are an 
outcome of a subtle interplay between chain fiexibility, 
size of the enclosure and steric interactions between DNA 
segments. It appears that the increase of the persistence 



length which would amount to the same thing as decreas- 
ing the size of the enclosure would promote the transition 
from folded toroid to a spool-like structure [54] . 

It thus seems that several salient features of the simula- 
tion results are in accord with our mesoscopic modeling 
of the confined polymer nematic ordering. First of all 
the robustness of the inverse spool configuration tran- 
spires from both types of approaches and points to its 
possible universality. Though ordered states discussed in 
this work never seem to point to the existence of folded 
spool configurations we do get other types of related or- 
dering with pronounced depleted regions and no polar 
symmetry. 

Though one gets the impression that molecular simu- 
lations of nematic polymer packing can incorporate more 
relevant details of molecular interactions and chain fiexi- 
bility, it should be stated quite clearly that the details of 
DNA-DNA interactions in particular are still not quanti- 
tatively known with great precision [21 . The mechanism 
of polyvalent-counterion mediated attraction keeps defy- 
ing proper understanding since apart from the charge 
of the counterions ^6\ it appears to be related also to 
ion-specific non-electrostatic interactions that are noto- 
riously difficult to parameterize and model [57 . The 
strong point of mesoscopic modeling as pursued in this 
work is that at least in principle one can include the 
experimentally measured density dependence of the free 
energy into the calculation directly by generalizing the 
quadratic term, Eq. (23), without invoking any particu- 



lar microscopic mechanism leading to the experimentally 
observed equation of state. This type of generalizations 
will be pursued in the future. 
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